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Abstract 

The increase of aqueous solubility of nonpolar compounds upon cooling and the cold denaturation 
of proteins are established experimental facts. Both phenomena have been hypothesized to be 
related to restructuring of the hydrogen bond network of water around small nonpolar solutes 
or hydrophobic amino acid side chains. However, an underlying physical mechanism has yet to 
be identified. We assume the solute particles and the monomers of a polymer interact via a 
hard sphere potential. We further assume that the solvent molecules interact via the two-scale 
spherically symmetric Jagla potential, which qualitatively reproduces the anomalies of water, such 
as expansion on cooling. We find that this model correctly predicts the increase in solubility of 
nonpolar compounds and the swelling of polymers on cooling. Our findings are consistent with the 
possibility that the presence of two length scales in the Jagla potential — a rigid hard core and a more 
flexible soft core — is responsible for both phenomena. At low temperatures, the solvent particles 
prefer to remain at the soft core distance, leaving enough space for small nonpolar solutes to 
enter the solvent thus increasing solubility. We support this hypothesized mechanism by molecular 
dynamic simulations. 



I. INTRODUCTION 



The increase of the solubility of nonpolar compounds in water on lowering the temper- 
ature, T, is a well known phenomenon which recently received attention from the 
molecular simulation community . The existence of closed demixing regions in the 
temperature-concentration phase diagrams characterized by the presence of a lower critical 
solution temperature (LCST) below which the components mix in all proportions is also well 
known for aqueous solutions of polymers and large organic molecules Band is explained b 



the difference in the thermal expansion coefficients of solvent and solute 
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The same effect, increase of solubility of the hydrophobic side chains on lowe ring T, may be 
relevant for understanding the remarkable phenomenon of cold denaturation, 
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a subject of considerable current theoretical 
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The thermodynamic framework of the solubility is well established [12j. The excess Gibbs 
free energy Ag of dissolving one mole of solute in x moles of solvent is defined as the difference 
between the Gibbs free energies of the solution and the sum of the Gibbs free energies of 
the pure solvent and the pure solute at constant pressure P and temperature T, 



Ag = Ah- TAs = Au + PAv - TAs, 



(1) 



where Ah, Au, Av, and As are the differences between the molar enthalpy molar potential 
energy, molar volume and molar entropy of the solution and pure substances. The main 
contribution to the entropic term As = Asq + As m ; x , is the entropy of mixing 



As mix ($) = -R 



ln$ + ^-lj ln(l-$) 



(2) 



where <& = 1/(1 + x) is the mole fraction of solute, R = N^ks is the universal gas constant, 
Na is the Avogadro number and ks is the Boltzmann constant. The term TAs m ; x dom- 
inates at high temperatures and thus favors mixing, but at low temperatures it becomes 
insignificant. However, the entropic term is not responsible for the solubility minimum as 
function of temperature. According to Le Chatelier's principle, if solvation is exothermic 
(Ah < 0) the solubility increases upon cooling at constant pressure. Conversely, if Ah > 0, 
the solubility increases upon heating. Thus the condition Ah = corresponds to the solu- 
bility extreme. Since d(Ag/T)/dT\p = —Ah/T 2 , where Ah = Au + PAV, the maximum in 
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Ag/T corresponds to the point at which Ah = 0, changing its sign from positive to negative 
upon cooling. 



These thermodynamic arguments were applied by Paschek et al., [171, [20j who showed 
that for the solution of Lennard- Jones polymer- chains in TIP5P water the enthalpic term 
Ah becomes increasingly negative at low temperatures and thus stabilizes the solution. 
The phenomenon of a solubility minimum has been discussed in terms of hydrogen bond 
restructuring in water at low temperature, which leads to the formation of water cages 



around small nonpolar molecules, monomers, or protein residues [2ll . |22| . |23| . Indeed, the 
PAv term becomes negative at low temperatures since, because the distances between water 
molecules increases upon cooling, there is more space for small solutes to penetrate the 
solvent. The Au term may also become negative, because introducing small nonpolar solutes 
inside water separates water molecules further apart and thus helps water molecules to form 
stronger hydrogen bonds. Without nonpolar particles, 5-coordinated molecules with weaker 
hydrogen bonds would be observed. 

In this paper, we attempt to elucidate the connection of the increase of solubility upon 
cooling with the anomalous thermodynamic properties of pure water. Recently it was shown 
that the spherically-symmetric Jagla ramp potential (Fig. [1]) explains anomalous thermody- 
namic and dynamic properties of water in terms of the presence of two length scales: a "hard 
core," corresponding to the first rigid tetrahedral shell of the nearest neighbors in water, and 
a wider "soft core," corresponding to the second (more flexible) shell of neighbors 24. l25|. 
We will see that the Jagla potential is capable of reproducing the solubility properties of 
water, and based on this result we will propose a mechanism for the solubility increase of 
nonpolar compounds in water on cooling. 

Using discrete molecular dynamic (DMD) simulations described in Section II, we show in 
Section III that the solubility of hard sphere solutes in the Jagla solvent has a characteristic 
minimum as a function of temperature similar to the minimum of the solubility of small 
nonpolar compounds in water. Hard spheres can be used to describe nonpolar compounds, 
since the van der Waals interactions among nonpolar compounds are an order of magni- 
tude weaker than the hydrogen bond interactions between water molecules and thus can be 
neglected in assessing the effect largely caused by the restructuring of hydrogen bonds. In 
Section IV we relate the solubility increase to the behavior of the excess thermodynamic 
functions. In Section V, we study the behavior of the beads-on-a-string homopolymer com- 
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posed of hard sphere monomers immersed in the same Jagla solvent. We show that the 
radius of gyration R g of the polymer has a sharp minimum as a function of temperature 
at low constant pressures, indicating that the system may display a LCST at low pressure. 
At higher pressure, the R g minimum becomes less pronounced, indicating the increase of 
solubility at intermediate temperatures and decrease of solubility at low temperatures. 



II. METHODS 



The interaction potential of the Jagla solvent particles U(r) is characterized by five pa- 
rameters: the hard core diameter a, the soft core diameter b, the range of attractive inter- 
actions c, the depth of the attractive ramp Ua, and the height of repulsive ramp Ur (Fig. [[]) 



261 ] . of which three are independent: b/a, c/a, and Ur/Ua- The most important of these 



parameters is the ratio of the soft core and hard core diameters, b/a, which must be selected 
close to t<2,/t\ % where t\ and r 2 are the positions of the first and the second peaks of the 



oxygen-oxygen radial distribution function for water. Following [25j, [26|, |27|, [28] we select 
b/a = 1.72, c/a = 3, Ur/Ua = 3.5. This choice of parameters produces a phase diagram 
which qualitatively resembles the water phase diagram with two stable critical points and a 
wide region of density anomaly bounded by the locus of the temperature of maximum density 
T M d [Fig. QJb)]. The slope of the coexistence line between the high density liquid (HDL) 
and the low density liquid (LDL) in this model is positive — unlike water, for which this 



slope is negative. In order to use the DMD algorithm [29|, [30J, |3JJ, we replace the repulsive 
and attractive ramps with 40 and 8 equal steps, respectively, as described in [25J (Fig. [[]). 
We measure length in units of a, time in units of a-sJm/U a, where m is the particle mass, 
density in units of a~ 3 , pressure in units of Ua/ci 3 , and temperature in units of Ua/^b- This 
realization of the Jagla model displays a liquid-gas critical point at T c % = 1.446, P c \ = 0.0417 
p c i = 0.102, and a liquid-liquid critical point at T c2 = 0.375, P C 2 = 0.243 p C 2 = 0.370 [25 ]. 
We model solute particles by hard spheres of diameter d , which we select to be equal to 
the hard core diameter of the Jagla solvent d = a. The dependence of the solubility on d 
is an important question and will be studied elsewhere. 
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III. CALCULATION OF THE HENRY CONSTANT 



We first study the effect of pressure and temperature on the solubility of hard-sphere 
solutes in the Jagla solvent. In order to do this, we create a system of N=1400 Jagla particles 
and 2800 hard spheres in a box L x xL y x L z with periodic boundaries. We fix L x = L y = 15a, 
and vary L z [Fig. QJc)] to maintain constant pressure and temperature using Berendsen 
thermostat and barostat. For T < T cl , the mixture of Jagla particles and hard spheres 
segregates and the Jagla particles form a slab of liquid crossing the system perpendicular 
to the ^-axis. For each pressure and temperature, we equilibrate the system for 500 time 
units and measure the mole fraction of hard spheres in the narrow slabs perpendicular to 
the z axis for another 500 time units. We find that this equilibration time is sufficient if the 
temperature for the two successive simulations is changed by less then 10%. In each slab we 
count the numbers Nj of Jagla solvent particles and N$ of hard sphere solute particles. 

We define a slab as a liquid slab if Nj > N and as a vapor slab if Nj < N , where N is a 
temperature and pressure dependent threshold, such that the liquid and vapor slabs form two 
continuous regions covering the entire system with the exception of a few slabs representing 
the boundary. To minimize the boundary effects, we exclude from the liquid phase the slabs 
whose distance to the closest vapor slab is less than 6a. The analogous criterion is applied 
to determine the vapor phase. Finally, we find the mole fraction $ = Ns/(Ns + Nj) of 
hard spheres in the liquid and vapor phases and compute the Henry constant fc#(T, P) = 
P$„(T, P)/$^(T, P), where Q v and $^ are the average mole fractions of hard spheres in the 
vapor and liquid phases respectively. 

Figure[2](a) shows the inverse Henry constant as a function of T for P = 0.1, 0.2, and 0.3. 
The behavior of the solubility follows the behavior of the inverse Henry constant [Fig. EJ^b)], 
since for T <C T C1 the partial vapor pressure of solvent is very low and, therefore, P v ~ P. 
Hence the solubility $£ « P/kn{T) is inversely proportional to the Henry constant and thus 
the solubility minimum as function of temperature practically coincides with the maximum 
of kn{T). One can see that Henry's law dkn{T, P)/8P\t = is approximately correct 
since kn(T,P) increases by less than 30% when the pressure increases by 200% from 0.1 
to 0.3. The most interesting feature of the behavior of the Henry constant is that it has a 
maximum as function of temperature at Tmh = 0.85 ± 0.10, indicating that the solubility 
of hard sphere solutes in a Jagla solvent has a minimum at this temperature and increases 
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upon cooling below Tmh- This behavior is similar to the behavior of the Henry constant of 
alkanes (C n H2„+2) in water, which increases upon cooling below T=373K [23} . Note that 
the T M h is much larger than the temperature of maximum density, T MD = 277K, for both 
water and for the Jagla model (for which Tmd ~ 0.5 at P = 0.1) 25]. 



IV. RELATION OF THE HENRY CONSTANT TO THERMODYNAMIC QUAN- 
TITIES 

To relate the behavior of the Henry constant to thermodynamic state functions, we must 
express the Henry constant in terms of Ag. At low solute mole fraction $^ = y^j (xe + ye) in 
the liquid phase, where xi and ye are the number of moles of solvent and solute respectively, 
the Gibbs potential G(xe,ye) of the solution can be approximated by the first two terms in 
the Taylor expansion of G(x£, ye) around ye = [0 

G e (x e , ye) = x e gi(0) + y t ^e - ^TAs mix ($ £ ), (3) 

where ge(0) is the molar Gibbs potential of pure liquid and fie is the chemical potential of 
solute in liquid. At low pressures, the analogous approximation can be made for the vapor 
phase 

G v (x v ,y v ) = x v g v (0) + y v fi v - y v TAs mix ($ v ). (4) 

Minimization of the sum of the Gibbs potentials of the two phases with respect to ye at 
constant y = ye + y v gives the equilibrium condition: 

exp — . (5) 



by 

§e V RT 

The chemical potential of solute in vapor is approximately equal to the molar Gibbs potential 
of the pure solute. Thus fie — fi v = Ag , where A g = Ag + TAs mix and Ag is given by 
Eq. (pQ). If we assume that the vapor is an ideal gas, then \i v = (J>o(T) + RTlnP, and thus 
Eq. ([5]) can be rewritten in the form of the Henry law: P v /&e — kn{T, P), where P v = P$„ 
is the partial pressure of solute in vapor and 

k H (T,P)=P eX p(jj[£) (6) 
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is the Henry constant which in the limit P — > does not depend on P. Thus the Henry 
constant has a maximum at the same temperature as Ago(T) /T. Since the entropy of mixing 
does not depend on temperature, this maximum coincides with the maximum of Ag(T)/T, 
which as we have seen above corresponds to the temperature at which Ah = 0. 

To study the behavior of each thermodynamic term in Eq. ([I]) at high and low pressures, 
we compare a mixture of Jagla particles and hard spheres with a given mole fraction $ and 
a given pressure P, with a pure solvent and pure solute at the same pressures. We simulate 
the system at constant P and $ in a cubic box with periodic boundaries with given numbers 
Ns and Nj of the solute and Jagla particles, respectively, such that Nj + N$ = 1728 and 
Ns/{Nj + N s ) = $. The systems were simulated at two pressures, P = Pi = 0.3 > Pq 2 
[path a of Fig. QT»] and P = P 2 = 0.1 < P C2 [path (3 of Fig. QTV)]. To find the entropy of 
the system as function of T and P, we change temperature from T = 10 to T = 0.3 with a 
small step and perform thermodynamic integration up to T = 10 using the first two terms 
of the virial expansion for T > 10 with our analytically computed values of the second virial 
coefficients for hard spheres and Jagla particles. We make sure that the system does not 
phase segregate at any intermediate temperature. 

Figure [3] shows the behavior of the dimensionless excess state function Au/RT, PAv/RT, 
Ah/RT, Aso/R, and Ag^/RT. One can see that Ah/RT becomes negative at the temper- 
ature Tmh(^)-P) of maximal Ag^/RT and hence, by Eq. ([6]), maximal Henry constant. We 
find T MH ($ = 0.2, P = 0.3) « 0.90 and T MH ($ = 0.1, P = 0.1) « 0.82 thus only weakly 
depending on $ and P. For both P values, Av becomes negative as the temperature drops 
below the liquid-gas critical temperature T Cl = 1.446. Since Am decreases with decreasing 
temperature in this region, Ah becomes negative at sufficiently low T. 

This behavior is consistent with the fact that the solute particles enter relatively large free 
spaces between the solvent particles which are kept apart not by the hard core but rather 
by the repulsive ramps, thus creating large negative Av [Figs, d^b) and [0(c)]. However at 
high temperatures, the entropic contribution to the free energy is important and the solvent 
particles in the solute are spread further apart than in the pure sovent. Thus the number of 
solute particles in the attractive range of the Jagla potential is reduced, leading to Am > 0. 
As the temperature decreases, the potential energy of the solution approaches its ground 
state, which is the same as that of the pure solvent, because the solute particles are small 
enough to remain in the area unavailable for the solvent particles at low temperatures due 
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to their repulsive ramp interactions [Fig. [U(b)]. 

At sufficiently low T, the excess enthalpy Ah of the solution becomes negative, indicating 
the increase of solubility upon further cooling. Close to T MD , the behavior at low and high 
pressure becomes different. For P = P 2 < P C 2 [Fig. E(b)], the pure solvent starts to expand 
upon cooling, thus increasing the free space for the solute particles, and further increasing 
the solubility. This is clearly indicated by the sharp decrease in Av. For P = P\ > 
[Fig. Eta)], there is no Tmd and the pure solvent collapses upon crossing the Widom line 
into an HDL-like liquid on the lower side of the Widom line. Thus Av starts to increase upon 
cooling, and eventually Av becomes positive. On the other hand Am keeps decreasing on 
decreasing T and becomes negative because the solute particles prevent the solvent particles 
from "climbing the repulsive ramps" upon compression and therefore help them remain near 
the minimum of the pair potential. 

To test this explanation directly, we compute for P = 0.3 and T = 0.5 the solvent-solvent 
density correlation function, gu (r), in the pure solvent and in the solution with solute mole 
fraction $ = 0.2, [Fig. 11(a)] . We see that the first peak in gn(r), corresponding to the 
solvent particles staying at the hard core distance in the pure solvent, significantly decreases 
in the solution. In contrast, the second peak, corresponding to the solvent particles staying 
at the soft core distance, increases. The behavior of the cross-correlation function guir) of 
the solute and solvent [Fig. E(c)] indicates that indeed the solute particles prefer to remain 
very close to the solvent particles in the region unavailable for solvent particles due to their 
repulsive ramp interactions. The behavior of the solute-solute correlation function g 22 {r) 
indicates a small decrease in the first peak and an increase in the second peak upon cooling 
[Fig. 11(d)] , corresponding to the increase of solubility. Overall, the structure of the solution 
is more pronounced than in the pure solvent case, indicated by negative Aso [Fig. [3(b)]. 

This behavior of the Jagla solvent with hard sphere solutes is analogous to the behavior 
of water molecules in the presence of nonpolar compounds. The hard core of the Jagla model 
corresponds to the rigid tetrahedron (first shell) of water molecules linked by hydrogen bonds 
to a central molecule. A particle at the hard core distance in the Jagla model corresponds 
to an extra (fifth) water molecule entering the first shell. This molecule may cause the 
formation of a bifurcated hydrogen bond [33j of higher potential energy than the normal 
hydrogen bond. The nonpolar solute molecules in water remain close to the first shell of 
water molecules, and thus prevent extra water molecules entering into the first shell and 
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forming bifurcated hydrogen bonds. Thus in water, solute particles help stabilize hydrogen 
bonds. The second shell of water becomes more structured by forming cages around the 
solute molecules without breaking any hydrogen bonds. 

For P = P 2 < P C2 , the Qxi{r) for pure solvent and solution are practically indistinguishable 
(not shown), with only slight decrease of the overall number of solvent particles within the 
attractive range in the solution compared to the pure solvent, corresponding to small Au > 0. 
The behaviors of #12 (0 and #22 (r) at low pressures (not shown) are practically the same 
as at high pressures. This means that for low pressures the increase of the solubility upon 
cooling is mainly because of the decrease of Av due to the small and even negative thermal 
expansion coefficient of the pure solvent. We expect that the relative strength of the two 
contributions, energetic and volumetric, to Ah strongly depends on the size of the solute 
particles. Accordingly, in water one might expect both mechanisms to be present, depending 
on P, and on the nature of the solute. 

To check the validity of our calculations of the Henry constant ku and the thermodynamic 
excess quantities, we compare fcn obtained directly by measuring the mole fraction of hard 
spheres in the slabs of Sec. Ill with kn computed using Eq. (El). Figure [2] shows that 
although the individual contributions of different terms to Ag behave in a complex way, the 
resulting behavior of the Henry constant is in good agreement with the direct simulations 
near the solubility minimum. Note the failure of thermodynamic predictions at high and low 
temperatures. The dramatic solubility drop predicted by Eq. ([6]) at P = 0.3 is because at 
fixed mole fraction, the amount of solute particles become insufficient to prevent the solute 
from collapse, and eventually Au also starts to increase upon cooling making Ah > at low 
temperatures. This predicted decrease of the solubility upon cooling is not observed in the 
direct simulations (Fig. [2]). The discrepancy between the thermodynamic predictions and 
actual behavior of the Henry constant is due to the fact that at low temperatures the mole 
fraction <£v > 0.3 is so high, that the linear Taylor expansion of Eq. ([3]) is no longer valid. 
Thus, in reality, there are still enough solute particles to prevent the solution from collapsing 
into an HDL-like liquid. In the next section, concerning polymer solution, we observe the 
decrease of solvent quality upon cooling at P > P C2 because the polymer monomers do not 
have sufficient local mole fraction to keep the solution from collapse. 
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V. POLYMER SWELLING UPON COOLING 



In order to relate the maximum of solubility of nonpolar compounds to the polymer LCST 
and protein denaturation, we next study the behavior of a polymer composed of M = 44 
monomers modeled by hard spheres of diameter do = a. We model covalent bonds by linking 
the hard spheres with the potential 



^bond(^) 



oo r < d\ 

d\ < r < d 2 , (7) 
oo r > d 2 



where the minimal extent of the bond is d\ = a and the maximal extent is d 2 = 1.2a. We 
simulate for t m time units the trajectory of the polymer at constant T and P in a cubic 
box containing iV = 1728 Jagla solvent particles with periodic boundaries. We calculate the 
polymer radius of gyration R g (tk), where 



^k) — ^ 1 — ^ '- (8) 

at equidistant times tk = kAt, where At = 1, and we calculate the rms radius of gyration 
R g , where R 2 g = Y%=i R 2 g(tk)/fn, and m = t m / At. 

Figure [5] shows the behavior of the radius of gyration of a polymer as a function of 
temperature for six values of pressure. Also shown for the P = 0.02 case is the inverse 
Henry constant. One can see that for small pressures, R g reaches a deep minimum at 
T = Tmr(-P), which occurs above the maximum of the Henry constant T = Tym(P). As 
pressure increases, the minimum becomes less pronounced and Tmr(-P) shifts to higher 
temperatures and eventually, at P = 0.4, almost disappears. According to the de Gennes- 



Flory-Huggins theory 34j], the radius of gyration is a function of polymer length and the 
solvent quality \: 

R g = £N 1 / 2 f(z), (9) 
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where z = (2% — l)iV 1//2 and 



\z 



l ' z z < 

1 z = 0, , (10) 

z 2v - x z > 



where z/ pa 0.588 is the Flory exponent for R g . Values of x < 1/2 correspond to a "poor" sol- 
vent, in which the polymer chain is collapsed and R g ~ V 1 / 3 . Values of x > 1/2 correspond 
to a "good" solvent where the polymer swells and R g ~ V 1 '. The value x = 1/2 corresponds 
to the "theta condition" with R g ~ V 1//2 . Thus for a polymer of a given length, the radius 
of gyration monotonically increases with solvent quality and can be used to calculate lines of 
the equal solvent quality in the P — T plane. We tested numerically that near the minimum 
of R g for P < 0.2, R g ~ V 1 / 3 , so the solvent quality is poor. For T < 0.5, R g is compa- 
rable with R gv = 4.77 in a vacuum, and hence the solvent quality is good. Thus the theta 
condition corresponds to some temperature between T = 0.5 and the value of T at which 
R g has a minimum. Accordingly, the LCST which differs from the theta point by a term of 
the order of l/y/M {jmJ is also located in this region. Assuming that R g = 3.5 for N = 44 
corresponds to the theta point, we construct a region in the P — T plane of approximate 
location of the LCST which merges with the upper critical solution temperature (UCST) at 
P = 0.12 and T = 1.25 (Fig. [6]). The region of the P — T plane below this line corresponds 
to the poor solvent in which sufficiently long polymer chains segregate from the sufficiently 
concentrated solution. This region lies below P C2 and significantly above the Tub line in the 
P — T plane. On the high temperature side, the region where R g < 3.5 is bounded by the 
liquid-gas critical point. 

Since the driving force of the folding of globular proteins is the formation of the hy- 
drophobic core, one can imagine that the regions in the P — T plane in which the various 
proteins can fold into their native states must have similar shapes. 

VI. DISCUSSION AND CONCLUSIONS 

Our work suggests a physical mechanism for the increase of solubility of nonpolar solvents 
upon cooling (Fig. [2]) connecting the cold denaturation of a polymer to the LCST. These 
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phenomena are related to thermodynamic anomalies of pure water because both thermody- 
namic and solubility properties are reproduced in the simple Jagla solvent model. Both of 
these are caused by the existence of two repulsive scales in the effective model potential, the 
hard core corresponding to the first tetrahedral shell of water molecules and the soft core 
corresponding to the second shell of water molecules. However, the temperature range in 
which cold denaturation and the solubility minimum occur does not directly coincide with 
the regions of other water anomalies, which are generally placed at lower temperatures. 

Interestingly, the region of poor solvent, corresponding to the sharp minimum in R g lies 
below the critical pressure of the second critical point. Above this pressure, the behavior of 
the solvent quality changes, being almost constant in the wide range of temperatures much 
above T c2 , but dramatically decreasing as we approach T c2 , indicating that the HDL in the 
Jagla model is a poor solvent, while the LDL is a good solvent. This behavior is in accord 
with the thermodynamic studies of Section IV where a sharp increase in Ag upon cooling 
near the Widom line is found. The physical reason for this is that LDL in the Jagla model 
has a larger free volume for the hard spheres to enter and is also less ordered than HDL, so 
entering of the hard spheres does not destroy this order. Whether this situation for water is 
the same as or the opposite is not clear, since LDL in water, although having larger volume, 
is more ordered than HDL. New simulations, underway, involving water models that have a 
liquid-liquid critical point in the accessible region and a chain of Lennard- Jones monomers 
may answer this question. 
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FIG. 1: (a) The discretized two-ramp Jagla potential captures much of the essential physics corre- 
sponding to the first and the second shells in water. The diameter of the hard core r\ = a and the 
diameter of the soft core T2 = b ~ 1.72 core determine the two essential length scales corresponding 
to the first and second shells of water, (b) Schematic P-T phase diagram of the two-ramp Jagla 
model of plot [IJa). Shown are the two isobaric paths simulated, Path a {P = P\ = 0.3 > Pc 2 ) 
and Path f3 (P = Pi = 0.1 < Pc 2 )- Also shown are the liquid-gas and liquid-liquid critical points 
(shown as solid circles), the corresponding coexistence lines (dashed lines) and the Widom line 
(solid line). Also shown is the locus of temperature of maximum density labeled Tmd (dotted 
line), (c) A snapshot of a simulation box for calculation of the Henry constant at Pi = 0.3 > Pq 2 
and T = 0.95 > Tc 2 - Black circles represent the hard cores of the Jagla solvent, gray circles 
represent hard sphere solutes. 
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FIG. 2: (a) Symbols indicate our simulation results of the inverse Henry constant of hard sphere 
solutes in the Jagla liquid solvent for three different pressures. Lines indicate our theoretical 
predictions based on calculations of excess thermodynamic state functions — see Eq. © and Fig. [3J 
Note the failure of the theoretical predictions for both high and low temperatures in the P = 0.3 
case caused by the fact that Eq. © is valid only in the limit of small liquid mole fraction $^ of 
solute in liquid, (b) The solubility of the hard sphere solutes in the Jagla solvent as function of 
temperature. The solubility minimum roughly coincides with the minimum of the inverse Henry 
constant. 
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FIG. 3: Excess thermodynamic state functions of hard sphere solutes in Jagla liquid solvent com- 
puted for (a,b) Pi = 0.3 > P C2 and $ = 0.2 and (c,d) for P 2 = 0.1 < P C2 and $ = 0.1. One 
can see that the temperature at which Ah = corresponds to the maximum of Ag/RT (arrow). 
Since PAv is only weakly T-dependent at high P, the decrease of Ah = Au + PAv upon cooling 
for P = 0.3 is mainly due to the decrease of the potential energy An, while at low pressure, the 
decrease of Av plays a major role. Note that as the system enters the region of negative thermal 
expansion coefficient of the pure solvent (T < Tmd)j for low P becomes exceedingly negative 
corresponding to the simple physical picture that as the solvent is cooled, the solvent particles 
"descend from the repulsive ramp" toward the potential minimum at r 2 = b so that the average 
distance between them increases, allowing more space for small solute particles to enter the sol- 
vent. At high pressure, both Ait and Av become positive at low T. This behavior is due to the 
rapid restructuring of the system as it enters the region of high density liquid above the positively 
sloped Widom line, which emanates from the critical point located at T c2 = 0.375 and P c2 = 0.243 
[28|]. This behavior leads to the rapid increase of Ago corresponding to the predicted decrease 
in solubility for P = 0.3 (Fig. [2]). However the calculated solubility in Fig. [2] does not decrease 
because, at this pressure, the actual solute mole fraction is so large that it is outside the region of 
applicability of Eq. © . 
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FIG. 4: Solvent-solvent correlation function gn(r) computed at P = 0.3 for pure Jagla solvent, and 
for the solution with mole fraction = 0.2 at (a) T = 0.5 and (b) T = 0.95. These plots support 
the hypothesis that in solution the solvent particles prefer to stay at the soft core distance = b 
and not to "climb the repulsive ramp" occupied by the solute particles, (c) The solvent-solute 
correlation function gi2 indicates that the solute particles prefer to stay at the hard core distance 
r± = a away from the solvent particles. For comparison we plot the correlation function of pure 
hard spheres at the same temperature T = 0.95 and pressure P = 0.3, showing a much smaller 
peak, (d) The solute-solute correlation functions 522 (r) at two different temperatures and the same 
pressure P = 0.3 show that at low temperatures the solute particles prefer to stay apart from each 
other, indicated by the decrease on cooling of the first peak and the increase of the second peak. 
This effect, that the solute particles clump less at low T, is in accord with the increase of solubility 
upon cooling (cf. Fig. 2). (e) Schematic of the Jagla solvent particles (black) surrounding a hard 
sphere nonpolar solute (white) of the same diameter a. At low T the soft cores (gray) of the solvent 
particles do not overlap. The distance between the solvent particles corresponds to the potential 
energy minimum r\ = b in Fig. [U[a) , thus the solvent particles form a cage around the hard sphere 
nonpolar solutes. (f) Schematic of the same situation but now for T ~ Tmh, the temperature of 
maximum Henry constant and hence minimum solubility. 
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FIG. 5: Dependence on temperature of the average radius of gyration of M = 44 polymer com- 
prising M = 44 hard sphere monomers for different pressures. The dashed vertical lines show 
Tmh(P) of the maximum of the Henry constant. We also show that for the lowest pressures the 
solubility minimum of the hard sphere solute particles roughly coincides with the R g minimum. 
We find the rms radius of gyration in a vacuum is R gv = 4.77. The evolution of the polymer chain 
in the solvent is much slower than in a vacuum, so to calculate a reliable R g , we must carefully 
estimate the error bars. In order to do so we divide the entire set of m measurements into m/n 
equal groups of n < m/2 subsequent values Rg(tk), and find the average for each group {R g )i, 
i = 1,2, m/n and the standard deviation a n of {R 2 g )i- If the values {R 2 g )i were independent, the 
error bars on their average could be determined as a n j \/m/n — 1, and this estimate would not 
depend on n. In fact, for small n, o n j ^Jvnjn — 1 is an increasing function of n, which reaches the 
plateau at n = r/At, where r is the correlation time. Accordingly, we determine the error bar as 
max„[cr T1 /y / m/n- 1], checking that this quantity indeed reaches a plateau for a given number of 
observations m. If this quantity does not reach the plateau, it means that the time of averaging is 
insufficient and we must increase m. 
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FIG. 6: (a) Contours of equal R g in the P-T plane, showing regions of good and poor solvent; the 
numbers denote the value of R g . The filled circles C\ and Ci indicate the liquid-gas and liquid- 
liquid critical points. Note that at low P, on decreasing T along path (3, one passes from a region 
of good solvent (swollen "denatured" polymers) to a region of poor solvent (collapsed polymers) 
and finally to a region of good solvent ("cold denaturation" ) . (b) Schematic illustration of the 
loci of the lower critical solution temperature (LCST) and the upper critical solution temperature 
(UCST) for a polymer chain such as studied here. The UCST delineates the region of high-T 
swelling ("heat denaturation of a protein") while LCST delineates the region of low-T swelling 
("cold denaturation of a protein") sampled by path (5. 
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